Hysteresis and re-entrant melting of a self-organized 
system of classical particles confined in a parabolic trap 
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A self-organized system composed of classical particles confined in a two-dimensional parabolic 
trap and interacting through a potential with a short-range attractive part and long-range repulsive 
part is studied as function of temperature. The influence of the competition between the short- 
range attractive part of the inter-particle potential and its long-range repulsive part on the melting 
temperature is studied. Different behaviors of the melting temperature are found depending on the 
screening length (k) and the strength (B) of the attractive part of the inter-particle potential. A re- 
entrant behavior and a thermal induced phase transition is observed in a small region of (k, i3)-space. 
A structural hysteresis effect is observed as a function of temperature and physically understood as 
due to the presence of a potential barrier between different configurations of the system. 

PACS numbers: 61.46.-W, 64.60.Cn, 75.60.Nt 
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I. INTRODUCTION 

The study of the properties of self-organized systems 
has increased dramatically in recent years. This inter- 
est originates from the possibility to control the forma- 
tion of patterns having an important impact on applica- 
tions that use large-scale self-assembly to create specific 
pattern morphologies. This kind of structures occur in 
systems from diverse areas including chemistry [l| and 
biology 2]. In physics, it was predicted that systems 
which fall into this morphological category are gener- 
ated by the competition between short-range attraction 
and long-range repulsion fs*]. This competitive interac- 
tion appears in many systems, such as magnetic mate- 
rials 0, [al , colloids [g, Q , and two-dimensional electron 
systems [a Q . Experimentally, colloidal systems are one 
of the most studied systems, which in combination with 
theoretical predictions, may lead to the design of novel 
soft materials and to an understanding of the glass and 
gel state of matter Recently, it was observed that a 
similar type of pattern formation can also arise in partic- 
ular classes of ultrasoft colloids with a strictly repulsive 
inter-particle potential [ll| . 

Besides presenting a rich variety of cluster types and 
showing an excellent model for technological applica- 
tions, colloidal systems have the added advantage of the 
facility to control the interaction between particles and 
of real time imaging of their configuration through video 
microscopy. A wide variety of studies of colloidal systems 
were performed in order to understand the structure and 
dynamics of different kind of systems, such as colloidal 
particles interacting through a short-range attractive and 
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long-range repulsive potential [13, Eli [3 • ^'^^ instance, 
a binary system of superparamagnetic colloidal particles 
that are confined by a two-dimensional (2D) water-air 
interface and exposed to an external magnetic field per- 
pendicular to the interface showed diverse stable config- 
urations [iBl- Moreover, the authors observed that clus- 
tering appeared only for one type of particles, instead of 
both types of particles. 

Recently, several models were developed having a small 
number of interacting particles in order to understand the 
behavior of colloidal systems as a function of temperature 
, [TtI . [isj . Different melting scenarios were studied ex- 
tensively in systems consisting of charged particles for 
a range of different inter-particle interaction and trap 
types, such as a system of binary charged particles con- 
fined by a circular hard- wall potential interacting by a re- 
pulsive dipole potential [l^ , confined by a parabolic trap 
potential and interacting by a Coulomb inter-particle po- 
tential [20l|, and non-confined particles with short-range 
attraction and long-range repulsion interaction [2l| . 

Motivated by the increased interest in the behavior 
of systems of particles that are characterized by a com- 
petition between short-range attraction and long-range 
repulsion, we analyze here the melting of a system com- 
posed of a finite number of classical particles interact- 
ing through a potential which is composed of a repulsive 
Coulomb and an attractive exponential term. The parti- 
cles move in a 2D plane that are confined by a parabolic 
trap. The zero temperature configurations of this sys- 
tem were studied by Nelissen et al. They observed 
several kinds of topological different configurations (e.g., 
ring and bubble configuration). But from this rise the 
question about the stability of those kinds of configu- 
ration against thermal fluctuation. This motivated us 
to study the melting of these ordered configurations and 
to analyze the effect of the interplay between the short- 
range and the long range interaction on the melting pro- 
cess. We found, that if one increases the temperature, 
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some of the bubble configurations exhibit a thermally 
induced structural phase transition and a remarkable re- 
entrant behavior. In addition, we found that this system 
exhibits hysteresis behavior for the mean radial displace- 
ment when we increase and decrease the temperature and 
the configuration goes through a thermally induced struc- 
tural phase transition. 

This paper is organized as follows. In Sec. II, we de- 
scribe the mathematical model and the numerical ap- 
proach to obtain the ground state configurations. In Sec. 
Ill we analyze the melting process for different inter- 
particle potentials. Temperature induced re-entrant and 
hysteresis effects are discussed in Sec. IV. Our conclu- 
sions are presented in Sec. V. 



II. NUMERICAL APPROACH 

We consider a 2D cluster with N classical particles 
interacting through a potential composed of a repul- 
sive Coulomb and an attractive exponential term as in 
Ref. |21), |23]. The particles are kept together by a 
parabolic potential centered at the origin. The general 
dimensionless Hamiltonian of the system is written as: 



H = 



N 



N 

E 



1 



(1) 



where = |ri| is the distance of the ith particles from 
the center of the parabolic confinement. The energy 
and the distance are in units of Eq — (mcjg rl/2) and 
ro — {2q^ /meujQy^'^ , respectively. Notice that the B and 
K parameters determine the exponential contribution in 
the hamiltonian, where B determines the strength and 
1/k the interaction range of the attractive part of the 
inter-particle potential (the third term in Eq. ([1])). Note 
also that the k parameter is inversely proportional to the 
range of the attractive part in the potential. Now the 
state of the system is determined hy B, k and the num- 
ber of particles N. Temperature is expressed in units of 
Tq = Eo/ks, where is the Boltzmann constant. 

The ground state configurations (T = 0) of the two- 
dimensional system were obtained by the Monte-Carlo 
(MC) simulation method (using the standard Metropo- 
lis algorithm 23|) extended with the Newton optimiza- 
tion method 2J]. The particles were allowed to reach 
a steady state configuration after 10^ simulation steps, 
starting from different initial random positions. In the 
same time, we calculate the frequencies of the normal 
modes of the system using the Householder diagonaliza- 
tion technique [l^l ■ The configuration was taken as final 
if all frequencies of the normal modes were positive and 
the energy did not decrease further. 

In order to understand the trajectory of each parti- 
cle correctly, we study the melting of the system using 
Molecular Dynamic simulation (MD). The temperature 
of system was increased from T = (ground state config- 
uration) with successive steps of ST and equilibrating at 



the new temperature during 10^ MD steps, (with a typ- 
ical step size of At = 0.001). After this equilibrium, the 
average energy was calculated, together with the mean 
squared radial displacement given by 
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where p is the average inter-particle distance at zero tem- 
perature. The symbol () stands for an average over typi- 
cally 10^ - 10^ MD steps after equilibration of the system. 

The melting temperature was determined through a 
Lindemann-like criterion, which has been widely used for 
2D finite size clusters. This criterion states that melt- 
ing occurs when (uj^) reaches 0.1 of the inter-particle 
distance at zero temperature [1^. But what is essen- 
tial is that melting is characterized by a rapid increase 
of the fiuctuation of particles when temperature reaches 
the melting temperature. 



III. MELTING 

In the following, we analyze the melting temperature 
and as an example we consider systems composed of 20 
and 30 particles and study its dependence on the B and k 
parameters. These two systems are the typical examples 
of the T = configurations consisting of rings or bubbles. 



A. Dependence on B 

In this section, we study the melting temperature, as 
a function of the strength of the attractive part of the 
inter-particle interaction (B), for TV = 20 and TV = 30 
particles and a fixed value of k. We observe from Fig. [1] 
that different B-regions exhibit a different melting tem- 
perature (Tm) dependence. For instance, in the case of 

= 20 and k — 6, we find that the melting temperature 
decreases when B increases (see Fig. lUa)), whereas the 
opposite behavior is found for large B values (see Fig. 
[IJb)). The behavior of Fig. [Ija) is in some sense a sur- 
prising result since B is the strength of the short-range 
attraction, and one may expect that the larger the at- 
traction strength, the more packed the particles are, and 
therefore the higher the melting temperature. This con- 
cept is correct since the attractive part of the potential 
is large enough to compete with the repulsive potential 
part, but on the contrary, the opposite effect happens be- 
cause the total inter-particle potential is purely repulsive 
in the case of the cluster of Fig. [Ija). In this situation, 
when the value of B increases, the repulsive potential 
decreases [see Eq. ([T])] in the same way as the melting 
temperature due to the decrease of the coupling among 
the particles. This behavior was also observed for the 
system with 30 particles and k — 4 [Fig. [^c)] . 

In Figs, mb) and[T](d), we show the dependence of the 
melting temperature as a function of B for the cluster 
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FIG. 1: The total melting temperature (Tm) as a function of B for a system with k = 6, A'^ = 20 particles (a) and (b), and a 
system with k = 4, N = 30 particles (c) and (d). We considered ring configurations (a) and (c), and bubble configurations (b) 
and (d). Symbols are the numerical results and the solid line is a guide for the eye. 



with N = 20, K — 6 and N — 30, k = 4, respectively 
in the large B-region. As seen, the attractive part of 
the potential is large enough to form small clusters, and 
therefore the melting temperature increases with B as 
a consequence of the increase of the attractive part of 
the inter-particle potential. Hence when the system is 
composed of bubbles, the melting temperature increases 
with B due to the increase of the coupling among par- 
ticles in the small bubbles. The behavior of Figs. [Ha), 
[TJc) and Figs. [Ub), [TJ^d) clearly shows the dependence 
of the melting temperature on B, which is opposite for 
pure ring configurations and bubble configurations. 



B. Dependence on k 

In this section, we present the melting temperature as a 
function of k, the range of the attractive part of the inter- 
particle potential, for = 20 and = 30 particles and 
a fixed value of B. It is important to remember that the 
K-parameter is inversely proportional to the range of the 
attractive part of the inter-particlc potential. Although 
the attractive range presented in Fig. H^a) is too small to 
agglomerate particles, it is significant to change the melt- 
ing temperature for small k- values. In other words, when 
the value of n is increased from k = 4.5 to higher values, 
the attraction between the particles decreases and conse- 



quently the Coulomb repulsion increases. Consequently, 
the particles become more packed and the melting tem- 
perature increases. For large k- values shown in Fig. [U^a), 
the attraction range is too small either to form small clus- 
ters or to influence the melting temperature. This behav- 
ior is due to the fact that the coupling among particles 
saturates for large K-values, that is, the particles cannot 
be more densely packed and therefore the value of the 
melting temperature becomes almost constant. 

For the situation when bubbles are present (see Fig. 
[21b)), we find that the melting temperature also increases 
with increasing k. This general behavior is due to a de- 
crease of the attractive part of the inter-particle potential 
when K increases. Moreover, we observe that for small 
K- values shown in Fig. Hfb), the range of the attraction is 
large enough to form bubbles but does not influence the 
melting temperature considerably. In other words, the 
particles are weakly coupled in small bubbles for small 
values of k and the small difference of k does not change 
strongly the coupling among particles and therefore the 
melting temperature. On the another hand, for large n- 
values in Fig. [2Ib), the different attraction range changes 
the melting temperature due to the increase of the cou- 
pling among the particles. In this situation, a small dif- 
ference of K increases the packing among particles and 
therefore changes the melting temperature. A similar 
trend in the melting temperature is found for other clus- 
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FIG. 2: The total melting temperature (Tm) as a function of k for a system with _B = 6, = 20 particles (a) and (b), and 
A'^ = 30 particles (c) and (d). We considered ring configurations (a) and (c), and bubble configurations (b) and (d). Symbols 
are the numerical results and the solid line is a guide for the eye. 



ters with 30 particles [Figs. [2l^c) and[2l^d) respectively]. 
We can see that the behavior of the melting temperature 
is in general the same for the different configurations, 
i.e., the higher the value of k, the higher the coupling 
among particles and consequently the higher the melt- 
ing temperature. However, the melting temperature has 
different regimes as a function of k for the same configu- 
ration. 



IV. STRUCTURAL BEHAVIOR 

Recently, colloidal systems exhibited several new and 
interesting features, such as re-entrant behavior [isl . 
[H, [2^ and a hysteresis effect [l^j. In this section, we 
show for a specific short-range interaction that our sys- 
tem can present a different re-entrant effect and a hys- 
teresis behavior as a function of temperature. 



In this section, we observed clearly the different de- 
pendence of the melting temperature with respect to the 
strength and the range of the attractive inter-particle po- 
tential. For different strength of the attractive potential, 
the behavior of the melting temperature is determined 
by the configuration of the system, i.e., pure ring config- 
urations or bubble configurations. On the another hand, 
the range of the attractive potential is important to de- 
termine the trend of the melting temperature for a spe- 
cific configuration. The bubble configurations presented 
in this section did not exhibit a structural transition from 
bubble to ring configuration when the temperature is in- 
creased. This is a consequence of the fact that the melt- 
ing temperature of the ring configurations are smaller 
than that of the bubble configurations, i.e., the bubble 
configurations are more stable than the ring configura- 
tions for the strength and the range of the attractive 
potential presented in Figs. [Hb) , [Hd) , [21Ib) and[2I^d). 



A. Re-entrant behavior 

The temperature dependence of the mean squared ra- 
dial displacement < u|j > is shown in Fig. EJa) for a 
system with B — 7 and k — 4. As we can observe, when 
the temperature increases from T = to T = 0.003, the 
value of < > increases considerably until it reaches 
a plateau and remains almost constant until T = 0.0053 
(dashed line). For T > 0.0053, we observe that the value 
of < > decreases rapidly before it increases sharply, 
indicating that the system melts. This re-entrant behav- 
ior was observed both in experimental ^, 26] and theoret- 
ical studies [II, 111] which, in the present case, is due to 
an increase of the stability of the whole system caused by 
an increase of the symmetry of the system configuration. 
In other words, the value of the mean radial displacement 
decreases after T = 0.0053 as a result of the change in 
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FIG. 3: (a) The mean radial displacement (< >), (b) the 
mean distance of each closest particle from the center of the 
confinement potential, (c) the mean radial displacement and 
(d) the angular intrashell displacement with respect to the 
center of mass of the small bubble for a system with TV = 20, 
5 = 7 and k = 4. 



tance of each particle in such a way that ri, r2, rs, r^v 
correspond to the first, second, third,..., TVth closest par- 
ticles from the center of the cluster, respectively. After 
that, we calculate the average of each closest distance 
< Ti > for each temperature and we present it as a func- 
tion of temperature in Fig. [3ljb). As can be seen, the sys- 
tem remains in the bubble configuration from T = to 
T = 0.0053. We will label this configuration as (2; 6(3)^) 
which means that there is one ring of 2 particles and 6 
bubbles of 3 particles [see the left inset of Fig. [Hb)]. For 
T > 0.0053, a particle leaves from the edge and goes to 
the center of the system. The configuration changes from 
the bubble (2; 6(3)-^) to the ring configuration (3; 17), 
which means that there is a ring of 3 particles and an- 
other with 17 particles. We observe from Figs. EJa) and 
^h) that the value of the mean radial displacement de- 
creases when the configuration changes from (2; 6(3)^) 
to (3; 17) and increases rapidly when the temperature 
approaches the melting temperature of the ring configu- 
ration. Therefore, this change of configuration is a struc- 
tural transition which stabilizes the system. This transi- 
tion occurs before the system is completely melted and 
is a thermally induced structural phase transition. This 
interest ing phenomenon was found in diverse previous 
studies [HjII^jIsO] but, in this case, it is the result of the 
increase of disorder in small bubbles when the tempera- 
ture approaches T — 0.0053. 

In order to better understand the local disorder of the 
particles, we present in Figs. [3^c) and ^d) the mean 
radial displacement and angular disorder of the particles 
of the small bubbles as a function of temperature. In 
particular, we calculate the radial and angular disorder 
of each particle with respect to the center of mass of the 
small bubble to which it belongs. Due to the rotation of 
the small bubble with respect to the confinement center, 
the center of mass of each small bubble is calculated in 
each MD step. The mean locally radial displacement is 
defined as 



'■R[B)I 



(3) 

where Nc and TV^c are the number of small bubbles and 
the number of particles in each small bubble, respectively. 
r|'" is the distance of the center of mass of the bubble 
from the center of the confinement potential and ps is 
the average distance between the particles of the same 
small bubble at zero temperature. 



the configuration of the system. 

In order to confirm that the re-entrant behavior is 
caused by the change in the configuration, we plot in Fig. 
[Hl^b) the position of the particles with respect to the cen- 
ter of the confinement potential as a function of temper- 
ature. Specifically, in each MD step we organize the dis- 



The angular disorder in the small clusters are stud- 
ied using the angular intrashell displacement calculated 
locally. Previously, this property was used to calculate 
the angular disorder in the whole system [l^, |23|, but 
in this case, we calculated the angular disorder in each 
small bubble, i.e., with respect to the center of mass of 
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FIG. 4: Phase diagram in {B, k) parameter space showing 
the re-entrant behavior for TV = 20 particles (solid symbols 
are the calculated values). The profile of the inter-particle 
potential for some relevant {B, k) values is shown in the inset. 



FIG. 5: The mean radial displacement for a cluster with = 
20, B — 7 and k = 4 for increasing (black ball and green 
triangle symbols) and decreasing temperature (red squared 
symbols) . 



the small agglomerate of particles virhich is defined as 

(4) 

where ii indicates the nearest neighbor in the same bub- 
ble and (pf) = 2tt/Nsc is the average density of particles 
in each small bubbles at zero temperature. 

As we can see in Fig. [3jc) and ^d), the mean ra- 
dial displacement and the angular disorder of particles 
increase simultaneously very dramatically in each bub- 
ble when the temperature approaches T = 0.0053 (dot- 
ted vertical line). This disorder permits a particle to 
overcome the potential barrier, between the center par- 
ticles and the edge particles, allowing to go to the center 
and changing the system configuration. Notice that the 
local disorder properties increase rapidly before the sys- 
tem changes its configuration from (2; 6(3)^) to (3; 17), 
that is, these properties change fast because of the in- 
crease of radial and angular disorder in the small bubble 
leaving initially the configuration unchanged. 

Re-entrant behavior was found in a small region of (k 
,i?)-space which is shown in Fig. U) We found that this 
interesting feature is the result of the inter-particle po- 
tential profile which is illustrated in the inset of Fig. 3] 
for some values of B and k. We observe that the inter- 
particle potential for the system with B = 7 and k = 4 
is not so repulsive as that for B — 6 and k — 5 and it 
is not so attractive as for B — 8 and k = 3. Therefore, 
systems which have values of B and k close to 7 and 4 re- 
spectively, are ideal to exhibit a re-entrant behavior due 
to the characteristic of agglomerating particles that are 
sufficiently weakly bound to allow particles to overcome 
the potential barrier when temperature increases, i.e., a 
thermally induced structural phase transition. 



The re-entrant phenomenon presented for the cluster 
with = 20 particles was also observed for other values 
of N, e.g. for = 30 particles. The general behavior of 
< w|. > shown in Fig. [Sja) is also observed for A^ = 30 
particles for a slightly different region in (k, i3)-space as 
well as the thermally induced structural phase transition. 



B. Hysteresis behavior 

Fig. [5] shows < > when we decrease temperature 
after the system has changed its configuration to (3; 17) 
for the same previous case, i.e., A^ = 20, B = 7 and 
K — A. Specifically, the mean radial displacement is cal- 
culated where the temperature is increased from T — 
to T = 0.008 (solid circles). In T = 0.008, the system 
does not melt but reaches a different configuration from 
the T = configuration. After that, we decreased tem- 
perature until T = (red square symbols) and increased 
it again up to T = 0.008 (green triangle symbols). Notice 
that at T = 0.008 the system reaches a different config- 
uration from the ground state one and that the < uj^ > 
behavior is very different when temperature is decreased 
indicating that the system gets stuck in a meta-stable 
state. This phenomenon is very interesting because the 
system exhibits a hysteresis effect in the mean radial dis- 
placement as a function of temperature when we decrease 
the temperature after a thermally induced phase transi- 
tion, i.e., after a change of configuration. This interest- 
ing behavior is a consequence of the fact that the (3; 17) 
configuration is stable, although it has a larger poten- 
tial energy than that of the (2; 6(3)'^) configuration, and 
is separated from the (2; 6(3)^) configuration by a high 
energy saddle point. 

In order to visualize this saddle point in the poten- 
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FIG. 6: The mean energy per particle as a function of the 
distance of the marked particle (particle in the inset indicated 
by a open circle) from the center of the confinement potential 
for the cluster with B = 7,it = 4, N = 20 particles and fixed 
temperature T = 0.0045. The configurations at the three 
local energy are shown in the inset as well as the trajectory 
of the marked particle (dashed line). 



tial energy landscape between the (3; 17) and (2; 6(3)-^) 
configuration, we fixed one particle from the edge of the 
(2; 6(3)^) configuration at T = 0.0045 and moved it in 
the direction of the center and at the same time allow- 
ing all the other particles to relax to their equilibrium 
positions. We calculated the average energy per particle 
of the system as a function of the position of the fixed 
particle [Fig. [5]. In this simulation we use Monte-Carlo 
(MC) technique in order to obtain the ground state con- 
figuration and to make an average of the particle energy 
of several MC steps (around 10^). That average energy 
per particle is presented in Fig. [S] as a function of the 
distance of the marked particle (open symbol in the in- 
set and this particle is moved along the thin dashed line) 
with respect to the center of confinement for a fixed tem- 
perature T = 0.0045. 

In Fig. [HI we observe that the average energy per parti- 
cle increases sharply when the marked particle is moved 
from r = 1.2 to r = 0.9. After that, the mean energy 
of the system decreases continuously until r ~ 0.3 and 
remains almost constant until r ^ 0.17. The (3; 17) 
configuration is found when r = 0.21. Notice that the 
(3; 17) configuration has the second lowest energy of the 
whole simulation and is stable. Thus, when we increase 
the temperature of the (2; 6(3)^) configuration such that 
there is enough thermal energy for a particle to overcome 
the potential energy saddle point, of T = 0.007, the sys- 
tem changes its configuration to (3; 17) and is locked into 



this lowest energy configuration. Consequently, the sys- 
tem remains in this meta-stable configuration even when 
temperature is decreased down to T = 0. 

V. CONCLUSION 

In this paper we investigated the dependence of the 
melting temperature of a system composed of classical 
particles, interacting through an inter-particle potential 
with a short-range attractive part and a long-range re- 
pulsive part, confined by a parabolic trap. The melt- 
ing temperature showed diverse behaviors as a function 
of the parameters that characterize the inter-particle in- 
teraction. In general, the melting temperature changes 
gradually as a function of the strength of the short-range 
interaction (_B), however, different behaviors of the melt- 
ing temperature were observed for pure ring configura- 
tions and bubble configurations as a function of B. The 
same trend of the melting temperature was observed for 
different range of the attractive part of the inter-particle 
potential (k) as well as for different configurations, that 
is, the melting temperature increases with an increase of 
the K parameter. However, different regimes of the melt- 
ing temperature appear in the same configuration as a 
result of a saturation of the coupling among particles. 

The mean radial displacement < u|j > showed a re- 
entrant behavior as a function of temperature. We found 
that this behavior is a consequence of a thermally induced 
structural phase transition which stabilizes the system 
before it melts. This structural transition occurs due 
to the rapid increase of the local disorder in the small 
bubbles of the system. This re-entrant behavior is found 
for a restricted set of values of the B and n parameter 
which define the attractive part of the potential. These 
values were shown in a phase diagram and showed that 
the re-entrant behavior is a characteristic of the inter- 
particle potential of the system. 

A hysteresis effect was observed in the structural and 
the dynamical behavior of the system as a function of 
temperature. It was shown that this behavior is a con- 
sequence of the existence of a high energy saddle point 
in the potential energy landscape between the two lowest 
energy configurations. 
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